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ABSTRACT 

We discuss the possibility of mapping interstellar clouds at unprecedentedly high spatial resolution 
by means of near-IR imaging of their scattered light. We calculate the scattering of the interstellar 
radiation field by a cloud model obtained from the simulation of a supersonic turbulent flow. Synthetic 
maps of scattered light are computed in the J, H and K bands and are found to allow an accurate 
estimate of column density, in the range of visual extinction between 1 and 20 magnitudes. We provide 
a formalism to convert the intensity of scattered light at these near-IR bands into a total gas column 
density. We also show that this new method of mapping interstellar clouds is within the capability 
of existing near-IR facilities, which can achieve a spatial resolution of up to ~ 0.1 arcsec. This opens 
new perspectives in the study of interstellar dust and gas structure on very small scales. The validity 
of the meth od has been recently dem onstrated by the extraordinary images of the Perseus region 
obtained bv lFoster fc Goodman I |)200!t1) . which motivated this investigation. 

Subject headings: ISM: clouds — ISM: structure — scattering — infrared: ISM — radiative transfer 



I. INTRODUCTION 

The spatial structure of interstellar clouds is the result 
of the interaction of supersonic turbulent flows with grav- 
itational and magnetic forces, which eventually leads to 
the formation of stars. Investigations of the cloud struc- 
ture can therefore provide constraints for models of star 
formation, as well as indirect evidence of the relative im- 
portance of turbulence, self-gravity and magnetic fields. 

Extended maps of interstellar clouds are usually ob- 
tained from i) the integrated intensity of emission lines 
of various molecular or atomic tracers, especially CO and 
HI, ii) the thermal emission of dust grains at far-IR and 
sub-mm wavelengths, and iii) the extinction of the near- 
IR light from background stars. Shortcomings of these 
methods are well known. A given molecular line, for 
example, is most sensitive to a limited range of gas den- 
sity, above the critical excitation density, and below a 
density where saturation sets in due to the large opti- 
cal depth. Furthermore, molecular abundances depend 
on complicated chemical networks affected by turbulent 
transport and depletion on dust grains. Thermal dust 
emission at far-IR and sub-mm wavelengths provides a 
more straightforward picture of column density, if the 
dust to gas ratio is constant and if the dust tempera- 
ture can be estimated. However, both dust temperature 
and optical properties of dust grains may have signifi- 
cant spatial variations, due to grain growth by coagula- 
tion and ice mantle deposition. Finally, reliable near-IR 
extinction maps can be computed only for regions at in- 
termediate Galactic latitudes and the uncertainty on the 
reddening is rather large as the spectral type of most 
background stars is usually unknown. 

In this Letter we discuss a new method of mapping 
interstellar clouds, based on the near-IR scattering of 
the general interstellar radiation field. This new method 
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avoids many shortcomings of the other methods men- 
tioned above. More importantly, it provides images of 
interstellar clouds at an unprecedented spatial resolu- 
tion, two orders of magnitude better than any single- 
dish line or continuum surveys. Near-IR scattering can 
be used to infer column densities within the approximate 
range of Ay = 1 to 20 mag. Even lower column den- 
sity can be probed with shorter wavelengths, using vi- 
sual and UV scattering. This study was motivated by 
th e extraordinary images of th e Perseus region obtained 
bvlFoster fc Goodman I 

Because of the very high spatial resolution (up to 
^ 0.1 arcsec with adaptive optics), near-IR scattering 
images offer exciting new perspectives in the investi- 
gation of the very small scale structure of interstellar 
clouds. 

2. TURBULENT CLOUD MODEL 

The density distribution of the model is the result of 
a three dimensional simulation of super-Alfvenic, com- 
pressible, magneto-hydrodynamic (MHD) turbulence, 
with rms sonic Mach number of the fiow, Ms = 10. 
The simulation is carried out on a staggered grid of 250"^ 
computational cells, with periodic boundary conditions. 
Turbulence is set up as an initial large scale random and 
solenoidal velocity field (generated in Fourier space with 
power only in the range of wavenumbers 1 < K < 2) 
and maintained with an external large scale random and 
solenoidal force, correlated at the largest scale turn-over 
time. The initial density and magnetic fields are uniform 
and the gas is assumed to be isothermal. Details about 
the numerical method are given in [Padoan & Nordlun^ 
(1999). 

The experiment is run for approximately 10 dynami- 
cal times in order to achieve a statistically relaxed state. 
The cloud model used in this work corresponds to the fi- 
nal snapshot of the simulation. The initial rms Alfvenic 
Mach number is A/a — 10. The volume-averaged mag- 
netic field strength is constant in time because of the im- 
posed flux conservation. The magnetic energy is instead 




Fig. 1. — Left: Projected density field from the super- Alfvenic MHD turbulence simulation with Mg 
in the central wavelength of the H band computed for the density field shown in the left panel. 



10. Right panel: Scattered light 



amplified. The initial value of the ratio of average mag- 
netic and dynamic pressures is {Pm)in/ {Pd)in — 0.005, 
so the run is initially super- Alfvenic. The value of the 
same ratio at later times is larger, due to the magnetic en- 
ergy amplification, but still significantly lower than unity, 
(^in)in/(-Pd)in = 0.12. The turbulcncc is therefore super- 
Alfvenic at all times. 

Supersonic and super-Alfvenic turbulence of an 
isothermal gas generates a density distribution with a 
very strong contrast of several orders of magnitude. It 
has been shown to provide a good description of the 
dynamics of molec ular clouds and of their highly frag- 
mented nature (e.g. 



Liar clouds and ot their highly h 
IPadoan et alHTOQl IMill Bool . 



3. SYNTHETIC MAPS OF NEAR-IR SCATTERING 

To compute the radiative transfer the data cube of the 
density distribution is scaled to physical units by fixing 
the length of the computational grid to L = 1.3 pc and 
the mean density to (nn) = 1000 cm^"^. However, the 
radiative transfer calculations apply to all values of L 
and (nn) satisfying the condition that their product is 
constant, L {nn) = 1300 pc cm~'^. When results are 
compared with observations, also the angular size of the 
models (or their distance) must be fixed. 

The model cloud is illuminated by an isotropic back- 
ground radiation, with intensities calculated according to 
iMathis. Mezerer fc Panagial l|1983D model of the interstel- 
lar radiation f ield (ISRF). Dust properties are taken from 
iDraind l)2003j) , and correspond to normal Milky Way dust 
with Ry = 3.1. We use the tabulated scattering phase 
functions that are available on the web^. 

The flux of scattered radiation is calculated with a 
Mont e Carlo program ijJuvela fc PadoaJ 120031 lJuvelal 
I2005j) . where sampling of scattered radiation is further 
improved with the 'peel-off' method (Yusef et al. 199j|). 
During each run, photons of one wavelength are simu- 

^ |http : //www ■ astro ■ pr inceton ■ edu/draine/dust/scat ■ html | 



lated and scattered intensity, including multiple scatter- 
ings, is registered toward one direction. The outcoming 
intensity is calculated separately for the central wave- 
lengths of the J, H, and K bands (1.25, 1.65, and 2.22 /zm) 
and for three directions perpendicular to the faces of the 
cubic model cloud. The result consists of maps of scat- 
tered intensity, where the pixel size corresponds to the 
cell size of the MHD simulation. The number of simu- 
lated photon packages is 4 x 10^ per wavelength. 

4. RESULTS 

The left panel of Figure shows a projection of the 
density field from the turbulence simulation described 
above. The corresponding image of scattered light in the 
H band is shown on the right panel of the same figure, 
for the model with mean gas density (nn) = 1000 cm~^. 
The H band scattering reproduces well the spatial struc- 
ture of the cloud model. The same is true for the J and K 
bands as well (not shown for lack of space). The scattered 
light tends to saturate only in the densest regions, at ap- 
proximately 5,10 and 20 magnitudes of visual extinction 
for the J, H and K bands respectively, as illustrated in 
the left panel of Figure [3 

Based on the radiative transfer along a single line of 
sight we may assume that the observed intensity / can 
be approximated with the function 
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-bN 
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where N is the column density. The coefficients a and 
b are positive constants defined for each band. They 
are related to dust properties and to the strength of the 
radiation field illuminating the cloud. To first order in 
N (small column densities), Eq. Ogives 

I = abN, (2) 

so the intensity is proportional to the column density. At 
large values of column density the relation is not linear, 
but can still be derived by comparing different bands. 
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Fig. 2. — Results from numerical modeling of light scattering. Left: J-, H-, and K-band surface brightness values as a function of A-y- No 
noise was added to these values. Right: Surface brightness of scattered light in the H band. Noise was added to bring the signal-to-noise 
ratio down to 10 (see text). 



For example, by writing Eq.^for the H and K bands, TV 
can be eliminated and the H band can be expressed as a 
function of the K band, 

7H=aHx(l-(l-:^)^). (3) 

flK 

The coefficients aj, an, and ok and the ratios 6j/6k and 
bn/bK can be determined by fitting this curve to the ob- 
servations (this is actually done in the three-dimensional 
(J,H,K) space). 

The b coefficients depend mainly on the properties of 
the dust grains. They should be well defined, since dust 
properties are believed to be very constant in the NIR. 
Their absolute values can be obtained from previous ob- 
servations of similar objects or by using stellar extinction 
measurements to fix the Ay scale in the current obser- 
vations. The a coefficients depend on the radiation field 
that illuminates the cloud. However, the spectrum of 
the incoming radiation can be directly estimated using 
observations at low extinction, and the field strength ap- 
pears basically as a multiplicative factor in the estimated 
column densities. 

The same deep observations carried out to detect the 
NIR scattering also provide color excesses for a large 
number of background stars. The intensity of the incom- 
ing radiation can be determined by comparing the stellar 
extinction with the amount of scattered light. All the 
coefficients can also be obtained from radiative transfer 
models based on the assumed dust properties and radia- 
tion field strength. Once the coefficients are determined, 
Eq. n defines, as a function of column density, a fixed 
curve in the three-dimensional (J,H,K) space. Column 
densities are estimated by finding the closest point on 
this curve to each observed (J,H,K) triplet and solving 
for N in Eq. [T] 

We have applied this method to the numerical cloud 
model described above. Noise is first generated so 
the signal-to-noise ratios for the mean intensities in 
the J, H,and K bands are 15, 10 and 7, respectively. 



The radiative transfer model gives aj^n.K = 1-04, 0.83, 
0.54x10-1^ ergcm-^sr-iRz-i and' 6j,h,k = 0.019, 
0.012, 0.0074 X 10"^° cm^, based on the true column den- 
sities and the scattered light observed toward one direc- 
tion of the cloud model. With these coefficients, the col- 
umn density for a different, orthogonal viewing direction 
is then estimated using the 'observed' surface brightness 
in that direction. Figure 13 shows the ratio between the 
estimated column densities and the true column density. 
The uncorrected values obtained assuming a linear rela- 
tion between the H-band intensity and the column den- 
sity are also shown. In regions with Ay ^10 mag the 
correction is almost a factor of three, or a factor ~2 if 
compared with uncorrected K-band intensities. Without 
correction, the estimated column densities are on the av- 
erage approximately 25% below the correct value. Fig.|21 
shows that the corrected estimates remain unbiased even 
at large extinctions. Ay ~ 10 mag. Their mean value 
deviates only 0.2% from the correct one, and the scatter 
around the mean is 10%. 

In Fig.|21i the scatter is real (no noise was added there) 
and due mostly to spatial variations in the field strength 
within the cloud. The field becomes more attenuated 
toward the cloud center, reducing the ratio of surface 
brightness to column density. In other words, the a co- 
efficients in Eq. ^ are functions of position. This is the 
dominant source of the scatter in Fig. 13 as well. Model 
calculations may be used to estimate the large scale vari- 
ations in the field strength, and to remove most of this 
scatter. A detailed radiative transfer calculation could 
also be carried out, where the cloud model is iteratively 
built to replicate the observations as close as possible. 
The estimated column densities would then only be sub- 
ject to uncertainties in the cloud structure along the line 
of sight. Other effects could be modeled as well, for ex- 
ample a non-isotropic radiation field. 

5. CONCLUSIONS 
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The average surface brightness of scattered light 
in the model with mean gas density (jt-h) — 
1000 cm-3 was 7.2 x IQ-i^, 4.9 x IQ-i^ and 2.8 x 
10~^^ ergcm"^ Hz^^ sr^^ in the J, H and K bands re- 
spectively. For a target signal-to-noise of 15, 10, and 7 
for the J, H, and K bands respectively, the necessary in- 
tegration times for the SOFT instrument on the 3.6 m 
NTT telescope become 33, 80, and 250 hours. How- 
ever, the instrument has 0.29" pixel size, and if reso- 
lution is degraded to 0.6" the quoted signal-to-noise ra- 
tios can be reached with less than 100 hours of integra- 
tion. With ISAAC/VLT, the same noise level per 0.15" 
pixel requires 22-1-72-1-196=290 hours. For a resolution 
decreased to 0.6", the integration time is reduced to ~17 
hours. 

At a distance of 300 pc, the angular resolution of 0.1" 
that can be achieved with adaptive optics, corresponds 
to a spatial resolution of ~ 30 AU. Imaging of near-IR 
scattering is therefore capable of probing the structure of 
nearby interstellar clouds at unprecedented spatial reso- 
lution, down to a physical scale of the order of the size 
of circumstellar disks, the gyroradius of GeV cosmic rays 
or the Kolmogorov turbulence dissipation scale. 

The right panel of Figure |3 shows the H band of our 
model. The image is composed of 250 x 250 pixels. After 



decreasing the resolution to 0.6", the same noise level and 
angular dynamical range is reached with ISAAC /VLT in 
~17 hours. This example shows that imaging of near-IR 
scattering is a promising method also for relatively large 
scale surveys of interstellar clouds, especially with the 
new generation of wide-field instruments, such as WF- 
CAM on UKIRT and WIRCam on CFHT. As discussed 
in the previous section, the method can provide a rather 
accurate estimate of the column density, preferentially at 
values of visual extinction between approximately 1 and 
20 magnitudes, and even below one magnitude if comple- 
mented with shorter wavelengths. It is therefore ideal for 
the study of low and medium density interstellar clouds, 
and less suitable for the interior of very dense protostellar 
cores. 
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